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ABSTRACT 

We propose a framework for understanding the fragmentation criterion for self- 
gravitating discs which, in contrast to studies that emphasise the ‘gravoturbulent’ 
nature of such discs, instead focuses on the properties of their quasi-regular spiral 
structures as derived from simulations. This interpretative framework is shown to be 
consistent with existing 2D and 3D numerical studies as well as with the 2D grid 
based and SPH simulations conducted here. We propose two evolutionary paths to 
fragmentation and argue that the correct simulation of each of these involves different 
numerical requirements: i) collapse on the free-fall time, which requires that the ratio 
of cooling time to dynamical time (/?) < 3 and ii) quasistatic collapse on the cooling 
time at a rate that is sufficiently fast that fragments are compact enough to withstand 
disruption when t hey encounter spiral features in the disc. We argue that the previous 
finding of Paardekooper et al. (2011) (in which 2D grid based simulations demon¬ 
strate a numerically converged fragmentation limit of (3 < 3 and which we reproduce 
here with both 2D grid based and 2D SPH simulations) is a consequence of the fact 
that such simulations smooth the gravitational force on a scale of H , the scale height 
of the disc. Such simulations thus only allow fragmentation via route i) above since 
they suppress the quasistatic contraction of fragments on scales < H ; the inability of 
fragments to contract to significantly smaller scales then renders them susceptible to 
disruption at the next spiral arm encounter. On the other hand, 3D simulations (along 
with 2D simulations that, with questionable realism, smooth gravity on smaller scales) 
indeed show fragmentation at higher /3 via route ii). We derive an analytic prediction 
for fragmentation by route ii) based on the requirement that fragments can contract 
sufficiently to withstand disruption by spiral features, basing this calculation on the 
properties of spiral structures derived from simulations. We find that this leads to a 
predicted maximum [3 for fragmentation of ~ 12 , in good agreement with all previous 
well resolved 3D simulations. We also discuss the necessary numerical requirements 
on both grid based and SPH codes if they are to model fragmentation via route ii). 


1 INTRODUCTION 


Observations have identified massive planets orbiting at 


large distances from their host stars (Marois et al. 20081. 


Such planets cannot have formed at their present location 
via the core-accretion mechanism and migration after for¬ 
mation at smaller radii is unlikely in some cases (Fabrycky 
& Murray-Clay 2010| ). An alternative explanation is that 
such planets formed near their present locations via direct 
gravitational collapse ( Boss| |2000). However, the exact disc 
properties necessary to form planets via the gravitational 
instability are still unknown. 

The gravitational stability of a disc has typically been 


understood in terms of the toomre Q parameter (Toomre 


1964), given by 


where c 3 is the speed of sound, E is the surface density of 
the disc and k is the epicyc.lic frequency (D for a Keplerian 
disc). Perturbation analysis of an axisymetric, infinitely thin 
disc indicates that Q > 1 is required for stability. 

When Q < 1 the gravitational instability creates spiral 
waves. The shock heating produced by these spiral waves 
stabilises the disc by increasing Q. This leads to a self- 
regulated Q ~ 1 state in which cooling is balanced by spiral 
wave heating. However, if the gas can cool efficiently, such 
a state is not established and over-densities instead collapse 
into bound objects. 

More than a decade of computer simulations have con¬ 
firmed that discs with Q < 1 are unstable, but can only 
fragment when cooling is efficient (Gammie|2001 Rice et al. 
|2005| |Meru fe Bate||2011| |2012| |Paardekooper et ah 


Q = 


ttGE 


(1) 


2003 


2011 ). 


Because the evolution of Q is driven most rapidly by 
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changes in temperature (Lod atop 0081, much effort has gone 
into defining the cooling rate necessary for discs to fragment. 
The simplest approach parametrises the cooling rate via a 
constant /3 as 


du 

dt 


u 

'W 


( 2 ) 


where u is the internal energy per unit mass. 

Early studies found that, for equations of state with 
(P and 7 = 5/3 , /3 > /3 cr it ~ 6 was sufficient to pre¬ 


vent fragmentati on (|Rice et ah 20051. It was later realised 
by Meru fe Bate| ( 201 1| ) that this result was not numerically 
converged, implying that fragmentation for slower cooling 
rates (higher (3) may be possible. Such non-convergence is 
unexpected, since even the lowest resolution simulations re¬ 
solve the disc scale height H, which is expected to be the 
most unstable wavelength. 

Several possible explanations of this non-convergence 


have been proposed. Lodato & Clarke (2011); Meru & Bate 


120121 suggested that unwanted artificial viscosity heating 


(which is resolution dependent) could stabilise discs against 


fragmentation. Rice et al. (2012) suggested that a modi¬ 


fication of the implementation of /? cooling in SPH could 


resolve the convergence issue. Paardekooper et al. (2011) 


found that when starting from commonly used smooth ini¬ 
tial conditions, a boundary layer develops between laminar 
and gravitationally structured flows and can induce frag¬ 
mentation. 

Fundamentally, the non-convergence of the fragmenta¬ 
tion boundary can only be explained if the results are cor¬ 
rupted by numerical effects, or if it is physically necessary 
to resolve length scale small than the Jeans’ length. This 
latter interpretation would be at odds with expectations 
based on the dispersion relation for gravitationally unsta¬ 


ble discs (Section 4.2I, which implies that the most unstable 
wavelength is ~ H and that the required amplitude for trig¬ 
gering collapse increases at smaller and larger spatial scales. 


Hopkins & Christiansen (20131 suggested that rare high am¬ 


plitude fluctuations on scale < H can lead to gravitational 
collapse and used this to argue for a stochastic fragmenta¬ 
tion model that could operate even at very large values of 
/?. Part of the motivation for the present study is thus to 
examine whether, in high resolution studies (with up to 30 
resolution elements across H), we see any evidence of this 
second mode. In Section [3] of this paper we present a range 
of new simulations examining the fragmentation boundary 
in 2D using both SPH and FARGO. We find that fragmen¬ 
tation is never initiated on scales <C H (in contrast to the 
hypothesis above), but that these scales become important 
in modelling the contraction of an unstable region to a size 
where it can survive disruption by spiral shocks. In Section 
[4] we expand on this argument, providing estimates based 
on characterisation of the quasi-regular spiral shock struc¬ 
ture in the simulations. Finally, in Section [5] we discuss how 
this relates to the results presented in the literature and its 
implications for the physical process of fragmentation. 


2 NUMERICAL METHODS 
2.1 Disc model 

We initialise all our simulations using a power law in surface 
density and temperature and a vertically isothermal Gaus¬ 
sian with scale height H. We fix the index of the surface 
density and choose the temperature scaling so that Q is ini¬ 
tially constant. That is, 


E = 


Mp 

nR? 


GM 

Ri 


-v p 
2 

R\ ~ p+3/2 qQ 0 T 
Ri 


(3) 

(4) 


where E is the surface density, M is the star’s mass, Mp 
the disc mass, Ri the inner radius of the disc and c s is the 
sound speed of the gas. q = Mp/M, Q o is the initial value 
of Q in the disc, p is a parameter to be specified and T is 
given by, 


r -i = f l«g(0 


(5) 


if p = 2 

\^(C 2 - p -l) if P 7^ 2 

where £ = R 0 /Ri and R a is the outer radius of the disc. The 
disc aspect ratio H/R is then given by 

2 -P 

( 6 ) 


H_ 

~R 


TQ 0 q I J? 
Ri 


while the ratio of the resolution scale h to the scale height 
H is, 


h 

H 


8-7T 


1 jd 


Ri) 


3(p —2) 


(7) 


^Nq 2 Q 2 0 TY 

where d is the number of dimensions. Note that the normali¬ 
sation of h/H will change by a factor of order unity between 
different types of code. 

We have slightly belaboured this point to emphasise the 
motivation for our parameter choices. The canonical choice 
in the literature is to choose p = 1. Equations [6] and [7] show 
that this choice results in a radial gradient in both aspect 
ratio and resolution, which breaks the otherwise scale free 
nature of the problem. To allow all radii in our simulations 
to be treated equally, we instead set p = 2. We also choose 
q = .2 and £ = 5 to reduce the computational expense of 


our simulations. Following Paardekooper et al. (20111 we use 


the equation of state P = KYP with the adiabatic exponent 
7 = 5/3. 

Our SPH simulations were performed using a mod¬ 


ified version of the popular code GADGET2 (Springel 
|2005[). The code was modified to include artificial conduc¬ 


tivity (Monaghan 19971, /? cooling, particle accretion and 


the correct treatment of softening with variable smooth¬ 
ing lengths (|Price & Monaghan 2007|). We also implemented 
the improved artificial viscosity method of |Cullen fe Dehnen| 
12010), which aims to minimise the amount of artificial vis¬ 


cosity present away from shock fronts, while still resolving 
shocks correctly. We set the strength of the artificial viscos¬ 
ity (and conductivity) to the values that produce the best re¬ 
sults in test problems where the correct result is known (e.g., 
shock tube test, Sedov blast wave, Kelvin-Helmholtz insta¬ 
bility). For artificial conductivity we use a con d = 1.0. For the 
standard implementation of artificial viscosity, the appropri¬ 


ate values are asPH = 1.0 and betasPH = 2.0 (Monaghan 
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1997). For the artificial viscosity method of Cullen & Dehnen 


12010), we use a max — 5.0, Ot-min — 1.0 and l = 0.05 ( |Cullenj 
fe Dehnen|2010 1. Note that although a ma x is quite high, in 


0.1 


practice the average per-particle value is only cxsph 
away from shocks. 

Grid based simulations were performed using the 
FARGO code (Masset 2000). FARGO uses a polar grid, 
which is logarithmically spaced in the radial direction and 
linearly spaced in the azimuthal direction. The numbers of 
radial and azimuthal bins were chosen so that each cell is 
approximately square (i.e. RA(f> « A R). 

All of the simulations in this paper are performed in 2D. 
This is partially to allow comparison with FARGO, which 
only operates in 2D with self-gravity, but also to maximise 
the simulation resolution ( h/H ), which scales as 7V -1 ' d ( d 
being the number of dimensions). 

After initialising each simulation as described above, 
we run each simulation for 700 td yn at the inner edge (10 
outer rotation periods) with /? = 30 to allow the disc to 
settle into the Q ~ 1 state without fragmenting. We then 
linearly decrease /3 from 30 to the desired value over the next 
700 tdyn- A simulation is deemed to be “non-fragmenting” if 
it does not fragment in a further 700tdj/n- This procedure 
follows |Paardekooper et al.| ( |2011[ > & jClarke et al.| ( |2007| ) 
and is designed to prevent fragmentation being induced by 
a boundary layer between regions of smooth and turbulent 
flow. 


2.2 Gravity in 2D 

The aim of two dimensional simulations is to capture the 
physics of the three dimensional system as accurately as pos¬ 
sible with a two dimensional representation. The treatment 
of the gravitational force requires particular care. In 3D, the 
gravitational potential is given by Poisson’s equation 

V 2 <I> = 4t rGp (8) 

In 2D, the code does not have access to the volume 
density and so we must obtain an expression in terms of E 
and other disc quantities. How best to do this depends on 
the details of the numerical code. FARGO, which calculates 
the gravitational force by directly solving Poisson’s equation 
uses the expression 

V 2 4> = 4ttGE5(2 - A) (9) 

where A is a softening factor typically set to some multiple 
of the disc scale height 77. If A = 0 then the code behaves 
as if A = h where h is the grid spacing. 

SPH calculates the gravitational force by summing the 
contribution from each particle (typically with the aid of 
a tree). To account for the vertical extent of the disc in 
2D, we calculate the gravitational force using a gravitational 
potential given by, 

N 

4? = —G m b (j>(\r — r b | , ft) (10) 

b= 1 

where mb is the mass of each particle h is some softening 


Fragmentation boundary using 2D SPH 

0.12 0.08 h/H 0.05 0.03 



Figure 1. SPH simulations as a function of cooling time ((3) and 
resolution (h/H), with circles denoting fragmenting and crosses 
non-fragmenting simulations, separated by the derived fragmen¬ 
tation boundary (solid line). The lower (blue) symbols are for 
simulations softened on scale H (the disc scale height), while up¬ 
per symbols are for those softened only on the resolution scale 

(h)- 


factor and 


fl//iG(|<?- fl 3 + §<? 4 ), 0 < g < 1; 

<t>=\ l/^(|g-3g 2 + |g 3 -lg 4 - I ^), 1 ^ q < 2; 

( 1/r 2 O 2 

(U) 

where q = r/ha (Price & Monaghan 20071. Adopting this 
form for (j> , it follows that 


V 2 <I> = 4t tG 


( — '■ 
\ 10 hn 


; 47rGp ( 


h si I 


( 12 ) 


where E is estimated using SPH interpolation with kernel 
length h,G ■ 

If we set h,G = hspH in SPH or A = 0 in FARGO, the 
right hand side of Equations [9] and [12] will increase with res¬ 
olution. That is, unless an explicit gravitational softening is 
provided in 2D simulations of self-gravitating discs, the grav¬ 
itational force will increase with resolution, suggesting that 
the value of the cooling time required to suppress fragmen¬ 
tation should increase with resolution. On the other hand, 
if the force is softened on the scale 77, one would expect 
simulations to converge because the suppression of gravita¬ 
tional effects on small scales means that these scales play no 
role in the fragmentation problem. To test this we ran all 
our simulations twice, once with no softening except on the 
resolution scale and once with softening on the scale 77. 


3 RESULTS 

3.1 SPH simulations 

To measure the fragmentation boundary, we performed a 
series of SPH simulations at different values of /3 using the 
gradually settled initial conditions described in Section [2] 
We repeated this experiment with the gravitational force 
softened according to /ig = 77 and /ig = hsPH■ All simula¬ 
tions were classified as “fragmenting” or “non-fragmenting” 
as described in Section [2] 

Figure [T] shows the results of these simulations with 
runs softened on 77 and h shown in blue and black respec¬ 
tively. Clearly the amount of gravitational softening used 
has a dramatic impact on the fragmentation boundary. The 
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Fragmentation boundary using FARGO 
cell size/# 

0.19 0.12 0.08 0.05 0.03 0.02 



Figure 2. FARGO simulations as a function of cooling time (3) 
and resolution (cell siz e/H), with circles denoting fragmenting 
and crosses non-fragmenting simulations, separated by the de¬ 
rived fragmentation boundary (solid line). The lower (blue) sym¬ 
bols are for simulations softened on scale H (the disc scale height), 
while upper symbols are for those softened only on the resolution 
scale. 


simulation softened on H show convergence to Bent ~ 3, 
which is consistent with the value found by earlier 2D stud¬ 
ies ( |Gammie|2001| ) & |Paardekooper et al.| ( |2011[ ). As well as 
fragmenting at much larger values of B , those simulations 
softened only on h show no signs of convergence. 


3.2 FARGO simulations 

To better understand the non-convergence of the SPH simu¬ 
lations and to verify that the SPH results are not an artefact 
of numerical technique, we ran the same simulations using 
the grid code FARGO. Once again, we softened the gravita¬ 
tional force on either the resolution scale of the simulation 
or on H. The simulations were initialised and classified into 
“fragmenting” and “non-fragmenting” in the same way as 
the SPH simulations. 

The results are shown in Figure[2] The runs with gravity 
softened on H are shown in blue while those without any 
extra softening are shown in black. The simulations which 
are softened on the scale of H are numerically converged to 
Bcrit ~ 3, in good agreement with previous 2D simulations 
and the results of our 2D SPH simulations. Once again, when 
we turn off gravitational softening the simulations do not 
converge and fragment at similar values of B to the SPH 
simulations at high resolution. 

For those simulations that fragment (in both SPH and 
FARGO), the fragments form from structures with length 
scales ~ H. This remains true even for the unsoftened sim¬ 
ulations that do not converge. That is, the non-convergence 
of our simulations without gravitational softening is not a 
consequence of instabilities triggered on scales significantly 
smaller than H , but of how resolution affects the evolution 
of regions of size ~ H once they become unstable. 


4 COLLAPSE CRITERIA 

The only difference between the converging simulations soft¬ 
ened on H and the non-converging simulations softened on 
h in Section [3] is on scales smaller than H. Although un¬ 
softened simulations are physically incorrect (as they fail to 


account for the disc’s 3D structure), it is at first sight sur¬ 
prising that processes on scales much less than the Jeans’ 
length (H for a Q = 1 disc) should affect the fragmentation 
process. 

For a disc to fragment, gravitationally unstable over¬ 
densities must be able to collapse into gravitationally bound 


objects that can resist disruption (Kratter & Murray-Clay 


2011) & (Paardekooper 2012). To determine the role of 


scales smaller than H in fragmentation and estimate when 
fragmentation is likely, we investigate the evolution of grav¬ 
itationally unstable over-densities in a Q ~ 1 disc using 
simplified models for the self-regulated disc structure. 


4.1 Free fall collapse 


The collapse of an over-density into a bound fragment can 
happen in two different ways. If the cooling rate is fast 
enough, any extra heat generated by compression of the col¬ 
lapsing over-density is radiated away. In this regime, the 
over-density can never become pressure supported and so 
collapses in free-fall. This will be the case whenever the cool¬ 
ing time is less than the free-fall collapse time. 

For a self-gravitating disc, the free-fall collapse time is 


given by Kratter & Murray-Clay (2011) 




1 


VGp 


- \j 2ttQ Irly 


(13) 


Therefore, an over density will undergo free-fall collapse 
whenever t coo i <tff, which is the case whenever 


B< \JYkQ m 3 


(14) 


In this regime, the collapse time tff -C t sp i (the time 
scale on which condensations encounter the next spiral 
shock: see below) for all reasonable values of Q. Disc frag¬ 
mentation is therefore inevitable and a self-regulated steady- 
state cannot form if over-densities collapse in free-fall. 


4.2 Pressure supported collapse 


If B > 3 then gravitationally unstable over-densities become 
pressure supported before they can completely collapse. Fur¬ 
ther collapse can only occur when the extra heat generated 
by collapse is removed, which happens on the cooling time 
scale. Such pressure supported clumps will only survive if 
they can collapse before being disrupted by their environ¬ 
ment. 

Although the Q ~ 1 quasi-equilibrium is constantly 
changing, it is still possible to define the average geomet¬ 
ric properties of the spiral waves that are created by the 
gravitational instability ( Cossins et al.|2009 |. In particular, 
we describe the spiral waves using an azimuthal wavenumber 
m and a radial wavenumber fc, which may vary as a function 
of R. 

Let us now consider a patch of disc that has become 
gravitationally unstable. On average, the longest it will sur¬ 
vive for before encountering a spiral shock is 


27T 

sp! m | fl — 12 


27T 

— T tdyn 

-Pl 


(15) 


where £ = |fi — H p | /fl and is the pattern speed of the 
spiral arms. 
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The core of our argument is that the survival of over¬ 
densities (and hence the ability of a self-gravitating disc 
to progress unstable over-densities to gravitationally bound 
fragments) depends on an over-density’s ability to collapse 
before being disrupted by encountering a spiral shock. To 
first order, collapse of an over-density takes place on the 
cooling time scale t coo i = fftdyn and over-densities live for 
at most tspi before encountering a spiral wave. Therefore, 
collapse and fragmentation requires that t coo i < t ap i, which 
is true if 


P< 


2tt 

mi 


(16) 


Provided that the tight-winding assumption holds 
(|fcJ?| m), it must be the case that 

m £ 


M = 


kH 


(17) 


where M is the Mach number of the spiral shock. 

Substituting Equation [IT] into Equation [16] implies that 

2 tt 


Bcrit = — = 

Pcrtt mi MkH 


(18) 


Given a dispersion relation, it is possible to re-write kH 
in terms of the Toomre parameter Q and the Mach number 
M. If we assume an infinitely thin, tightly would disc, the 
standard dispersion relation can be re-written as (Bi nney fc] 
Tremaine|2008 1 , 


2/-2 u 2 j 2 2H\k\ 
mi = H k -—-b 1 

substituting in Equation [IT] and solving for kH gives 


kH = 


-1 + y/l + Q 2 (M 2 -l) _ Q 
Q{M 2 — 1) ~ 2 


(19) 


( 20 ) 


ffcrit — 


( 21 ) 


where the last equality holds when M — 1<1. 

Substituting Equation |20| into Equation 1 1 8| gives 

47T 

MQ 

where Q is to be determined numerically. 

Although it may seem that Q is unconstrained in this 
model, it is actually equivalent to setting the most unstable 
wavelength. For example, if the most unstable wavelength 
is similar to that of the axisymmetric disc (where kH = 
1 /Q), then it must be the case that Q = \/2. Numerical 
simulations typically find that 1 < Q < V2- Together with 
the expectation that Mm 1 , this leads us to predict that 
/3 C rit ~ 9 — 13 for well resolved simulations. 

Equations [18]&[21] are important as they provide a sim¬ 
ple analytic estimate of ff cr it (the first, to our knowledge). 
That said, the arguments made in reaching this estimate are 
all “first-order” and so we expect Equation[2l]to be accurate 
only up to a factor of order unity. Furthermore, Equation 
|19| assumes an infinitely thin disc. If we use a dispersion 
relation that accounts for disc thickness (see Equation 31 
of Cossins et al. (2009)), the resulting estimate of ff cr it is 


reduced slightly. More explicitly, 

Pcrit = MQ {1 ~ f ] (22) 

and the most unstable wavelength changes so that we expect 

Q = 1. 


4.3 Resolution requirements 


To understand the resolution requirements of resolving the 
collapse process described in the previous section, we need 
to know how “collapsed” a clump needs to be in order to 
survive a collision with a spiral arm. The average amount 
of energy per unit mass lost to cooling between spiral wave 
encounters is roughly ut ap i/t CO oi- If shock heating balances 
cooling (as it must in the Q ~ 1 state), material will gain 

2 t r 

An — ut s pi/t CO oi — (23) 

miff 


in internal energy per encounter with a spiral. 

In order for a clump to survive a collision with a spiral 
arm, it must still be smaller than its initial size, once this 
extra A u energy has been deposited. That is, collapse begins 
with a patch of size ~ H, which then contracts until it is 
hit by a spiral wave. This spiral wave deposits energy in the 
clump, causing it to expand. If the expanded size of the post¬ 
shock clump is greater than the size when collapse began, 
the clump will no longer be unstable and will not survive. 

The binding energy of the clump with size x, per unit 
mass, is given by, 


_ 2 GM C 

G-bind — 

6X 


(24) 


and so the post-shock clump size, X 2 , is given by, 


2 GM C _ 2 GM C 2t tu 

3x2 3x miff 

If we require that X 2 < H, this implies that 


(25) 


X 

H< 


1 + 


37 tQ 


mi'ri'r -1)/3 


(26) 


For reasonable values of 7 = 5/3., Q = 1 and ff = 6 , 
this yields the requirement that x/H < .25. Furthermore, 


if we derive a simple expression for x(t) (see Kratter & 
|MurrayIC lay] ( [2011 ) for example), substitute x{t ap i) for x 
using Equation 15 and solve for the ff (for the same values 
of 7 and Q), we find that fragmentation occurs when 


ff< 


2.2tt 

mi 


(27) 


a requirement that is remarkably similar to Equation |18[ 
which we obtained from simple time scale arguments. The 
fact that Equations [27] and [18] are so similar means that 
requiring t CO oi < t sp i ensures that the clump has contracted 
sufficiently to satisfy Equation [25] and hence survive. 

However, the derivation of Equation [26] has required a 
number of simplifying assumptions regarding the non-linear 
evolution of the unstable over-density. The true utility of 
Equation [26] is not in providing a more precise constraint on 
ffcrit than Equation [l 8 ] (both are subject to an uncertainty 
of order a factor of 2 ), but in understanding what physical 
scales and processes need to be resolved in order to accu¬ 
rately model the non-linear evolution of the clump and the 
formation of fragments. 

More specifically, we find that to survive a collision with 
a spiral wave, a fragment must have x/H <g 1. The exact 
value of x/H required depends on ff and like all other ar¬ 
guments in this section is uncertain to within a factor of 
order unity. Nevertheless, this reasoning shows why resolv¬ 
ing scales smaller than H is necessary to accurately model 
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the physics of disc fragmentation. It is necessary to resolve 
the physical scale of the initial instability and the scale to 
which that instability must shrink in order to survive colli¬ 
sions with spiral waves. 


5 DISCUSSION 


Most simulations that have investigated the fragmenta¬ 
tion boundary in 2D have found a numerically converged 
fragmentation boundary of P cr u ~ 3 ( jGammie 2001 

Paardekooper et al. 2011). These 2D simulations softened 
their gravitational force on ~ similar to our simulations 
in Section [3] which also found P cr u ~ 3. The only major 
exception to this is Meru & Bate (2012), who used a very 
weak gravitational softening of 3 x 10 ~*H and obtained a 
much higher fragmentation boundary and much weaker nu¬ 
merical convergence. The fragmentation boundary also ap¬ 
pears to be independent of 7 in softened 2D simulations, 
with simulations using 7 = 2 (Gam mie 2001 ) and 7 = 5/3 
(Paardekooper et al. (2011 1 , Section|3j) all yielding p cr it ~ 3. 

In contrast to this, all 3D simulations performed to date 
fragment at well above /3 = 3 ([Rice et al.||2005l [M2| [Mi) 


Cossins et al.|2009 Meru fe Bate| 20 lT 20121 . The fragmen¬ 

tation boundary in these simulations also shows differing 
degrees of numerical convergence ( jMeru fe Bate|2011[ |2012| 
Rice et ah||2014 l and have been convincingly shown to de¬ 
pend on 7 ( Rice et al.|2005 |. 

A consequence of softening gravity on the scale H is that 
the force between mass elements separated by less than H 
goes to zero (to first order). This suppression of the grav¬ 
itational force on small scales means that any clump that 
becomes pressure supported on a scale l < H lacks the in¬ 
wards gravitational pull necessary to drive it to contract fur¬ 
ther. I 11 effect, softening of small scale gravitational forces 
causes pressure-supported collapse to stall. Because of this, 
fragmentation via pressure-supported clumps is suppressed 
in 2D simulations with gravitational softening. 

It is only when t coo i < tff that pressure is never able 
to impede collapse and simulations softened on the scale H 
can fragment. In this regime, softening on H cannot impede 
collapse since the clumps enter the softened regime in free 
fall and so it is not necessary that they experience continued 
inward acceleration to drive collapse. Since tff~ 3td yn this 
occurs whenever P < 3, in excellent agreement with the con¬ 
verged fragmentation boundary identified in our simulations 
and in the simulations of Gammie ( 2001 ) & Paardekooper 


et al. (2011). If the amount of gravitational softening is re¬ 


duced, the fragmentation boundary will increase as the scale 
at which pressure supported collapse will be stalled will de¬ 
crease (M eru fe Bate|2012 l. 

If gravitational softening is set by the resolution scale, 
then the gravitational force will be given by, 


V 2 $ = 4ttG(E//i) 


(28) 


1 |Gammie| (2001 ) does not use an explicit gravitational soften¬ 
ing, but is forced to omit small wavelength modes in his calcu¬ 
lation of the gravitational force (which uses a Fourier transform 
of Poission’s equation). Discarding these modes is equivalent to a 
gravitational softening on the scale ~ 0.3 H. 


where h ~ IV -1 / 2 . That is, the gravitational force will con¬ 
tinue to increase with resolution. This increased gravita¬ 
tional pull will make fragmentation easier as the resolution 
increases, as is demonstrated in Figures [l] and [2] Even if 
a converged value of the fragmentation boundary can be 
reached with gravity softened on the resolution scale, its 
physical interpretation is unclear (Muller et al.|20l2l. 


Rice et al. (20051 showed that the fragmentation bound¬ 
ary varied as ( 7(7 — 1 )) _1 in 3D simulations. However, the 


softened 2D simulations of Gammie (20011, Paardekooper 
et al. ( 201 1 |) and those performed here all agree on the 


fragmentation boundary, despite using significantly differ¬ 
ent values of 7 . This finding is also consistent with different 
types of simulations (softened 2D, realistic 3D), probing dif¬ 
ferent modes of fragmentation (free-fall collapse, pressure 
supported collapse). Softened 2D simulations are only able 
to probe fragmentation by free-fall collapse, which does not 
directly depend on 7 (since the free-fall time is independent 
of 7 ). 3D simulations also permit fragmentation via pressure 
supported collapse, which does depend on 7 (see Equation 
26|. 

It is possible that a more sophisticated model for ap¬ 
proximating the gravity in two dimensions, such as the tech¬ 
nique proposed by [Muller et al. (2012|, would reduce the 
fragmentation suppressing effects of gravitational softening. 
However, it is inevitable that some differences with the full 
3D solution to Poisson’s equation will remain. In particu¬ 
lar, models that assume a vertical structure for the disc are 
likely to be least accurate at modelling the quasi-spherical 
over-densities that must be modelled for fragmentation. As 
some gravitational softening is necessary to approximate the 
three dimensional force in two dimensions, any two dimen¬ 
sional simulation must either fail to model the gravitational 
force 100 % correctly or suppress fragmentation via the col¬ 
lapse of pressure supported clumps. 

Three dimensional simulations do not suffer from this 
limitation and so are able to model both pressure supported 
and free-fall collapse. However, there remains some disagree¬ 
ment as to the converged value of the fragmentation bound¬ 


ary in 3D (Meru & Bate 2011 2012 Paardekooper et al. 


[MI] |Lodato fe Clarke|| 2011 [ |Rice et al.||2012[|2014| ). The 
model for pressure supported collapse developed in Section 
[4] estimates /3 cr it ~ 12. Given that we do not expect this 
estimate to be more accurate than a factor of 2 , this does 
not allow us to discriminate between the /3 cr it ~ 8 favoured 


by Rice et al. (2014 1 and the p cr u = 20 — 30 of Meru & Bate 


12012). However, we argue in Section 5.2 that the quasi¬ 


regular nature of the spiral structure in self-gravitating discs 
lead us to expect that fragmentation should not be possi¬ 
ble at value of p significantly higher than the latter range 
(Pcrit = 20 — 30), even allowing for the possibility of stochas¬ 
tic fragmentation. 


5.1 Numerical technique specific issues 

A fixed grid, such as the one used by FARGO, is only able 
to follow a contracting clump down to the resolution scale 
of the simulation. Given this, care must be taken that the 
resolution scale is much smaller than the smallest scale of 
physical importance. We have shown in Section [4] that this 
length scale can be significantly smaller than the disc scale 
height H. Failure to resolve small enough scales will suppress 
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Dependence of fragmentation in self-gravitating accretion discs on small scale structure 7 


fragmentation by preventing collapse from proceeding to the 
point where it can resist disruption by interaction with spiral 


By contrast, the adaptive nature of SPH ensures that 
so long as the resolution is high enough to resolve the initial 
instability ( h/H < 1), it will continue to resolve the clump 
as it contracts. However, modelling the collapse of pressure 
supported clumps is challenging in SPH because of the effect 
of particle noise. Any random “thermal” motions acquired 
by SPH particles will supply an additional pressure that will 
prevent collapse. Unlike pressure derived from internal en¬ 
ergy, this pressure due to particle noise can not be removed 
by cooling. 

Failure to provide adequate artificial viscosity at shock 
fronts is a well known source of particle noise in SPH. Strong 
shear, as is present in differentially rotating discs, can also 
lead to particle noise. Consistent with this interpretation, 


Meru & Bate (2012) found that lowering the strength of the 


artificial viscosity creates particle noise and inhibits frag¬ 
mentation. It is also worth noting that there is abundant 
evidence from tests with analytic solutions that the sup¬ 
pression of particle noise at shocks also requires a rather 


large value of oisph (> 0.7) Lattanzio et al. (19861; Lorn 
|bardi et al.| ( |1999] ); |Thacker et al.| ( |2000| ) . This is significantly 
larger than the value of osph = 0.1 this is commonly em 


ployed in disc fragmentation calculations l 

Rice et al. 

2005 

Cossins et al. 12009 Meru & Bate| 2011; |Forgan et al. 

2011 

Rice et al.|20141. Unfortunately, such a high value of a sph 


produces significant heating away from shock fronts unless 
the resolution is very high and this can suppress fragmen¬ 
tation, as has been discussed extensively in the literature 
(Meru & Bate|2011 2012 Lodato & Clarke|2011 Rice et al. 
|2014[ ) 

Therefore, special care must be taken to ensure that 
shock-fronts are well modelled and particle noise suppressed 
in order to accurately model fragmentation. Unfortunately, 
doing so greatly increases the number of particles required 
to prevent excessive artificial viscosity heating away from 
shock-fronts (although this computational cost can be re¬ 
duced by using modern SPH techniques such as improved ar¬ 


tificial viscosity triggers (Cullen & Dehnen 20101 and kernels 
with increased neighbour number (|Dehnen fe Aly|2012[|). 


tions in spiral geometry, it can survive longer than would be 
expected from the mean geometry and collapse enough to 
fragment. 

Precise numerical studies are required to quantify the 
exact likelihood that this random process will lead to frag¬ 
mentation, which is beyond the scope of the present study. 


5.3 General implications 


This paper has argued that two dimensional simulations can¬ 
not be trusted to model the fragmentation process correctly. 
Nonetheless, we are now in a position to provide some im¬ 
portant upper limits on when fragmentation can occur. 

We argue that fragmentation is impossible if the initial 
instability cannot contract significantly before it encounters 
a spiral arm. Since collapse occurs on the cooling time scale, 
this occurs whenever 


271 - 

tcool ^ T^dyn 

mt. 


(29) 


Furthermore, ml ; can be related to Q and the Mach number 
via the dispersion relation to give, 


47T 

tcool MQ ^ dyrL 


12 t d y 


(30) 


Equation [30] is uncertain to within a factor of a few. 
Nevertheless, we can confidently say that our model predicts 
that fragmentation will not occur whenever /3 is greater than 
a few tens. 

There is some possibility that stochastic fragmentation 
can lead to fragmentation at higher cooling times. Over the 
self-gravitating lifetime of a disc, this may shift the “effec¬ 
tive fragmentation boundary” a factor of two or so higher 
(although detailed numerical studies are needed to confirm 
this). 

The opacity regime in the outer parts of protoplanetary 
discs leads to /3 oc R ~ 9 ^ 2 ( Clarke|2009 Cossins et al.|2010 


|Paard ekoopcr 2012). Since /3 is such a strong function of 
radius, a factor of a few uncertainty in /3 cr it makes little 
practical difference to planet formation theory. That is, the 
upper bounds provided by Equation [30] already place the 
allowable radii for fragmentation with sufficient precision. 


5.2 Stochastic fragmentation 

It was shown by |Paard eko oper| ([2012]) that even when sim¬ 
ulations settle into the Q ~ 1 state without fragmenting 
immediately, they can still fragment stochastically over the 
lifetime of the disc. |Paard e kooper| (i2012| argued that this 
occurred when a clump “got lucky” and survived for long 
enough to resist disruption and form a fragment. This pro¬ 
cess can also be understood within the framework of quasi¬ 
regular, but intermittent, spiral structure. 

The expression which relates the geometry of the spiral 
waves to the fragmentation boundary, m£, describes the av¬ 
erage geometry of the disc. The spiral arms present at any 
particular time will fluctuate about this average geometry. 
When /3 > pcrit, the average time between spiral wave en¬ 
counters is shorter than the collapse time of an over-density. 
As such, the disc is resistant to fragmentation on average. 
However, if the gravitational instability happens to be trig¬ 
gered in a part of the disc with the right random fluctua- 


6 CONCLUSIONS 

In this paper we have presented a new interpretation of the 
fragmentation process in discs, which emphasises the role of 
the quasi-regular spiral structure. In this interpretation, a 
disc can only fragment when an over-density is able to con¬ 
tract to a size small enough to survive being disrupted when 
it encounters a spiral arm. We show that this requirement 
means that simulations of fragmenting discs need to resolve 
both the initial instability at scale H and the subsequent 
collapse of the resulting over-densities to a size x <C H. 

We have further shown that the fragmentation bound¬ 
ary obtained in 2D simulations depends strongly on the type 
of gravitational softening employed. We have found that al¬ 
though gravitational softening on the scale H is necessary 
to account for the vertical structure of the disc and obtain a 
converged fragmentation boundary, it also suppresses frag¬ 
mentation by preventing the collapse of pressure supported 
clumps. As such, it is perhaps unlikely that any 2D model 
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of a self-gravitating disc will be able to capture all the pro¬ 
cesses necessary to model disc fragmentation in its entirety. 

The necessity of modelling the contraction of a clump 
from instability through to disruption resistant fragment 
also leads to additional numerical requirements for simula¬ 
tions of fragmenting discs. For codes employing a fixed-grid, 
it imposes the requirement that scales x <C H be resolved 
so that clumps can collapse enough to resist disruption by 
spiral waves. For SPH calculations, where this resolution re¬ 
quirement is automatically satisfied if H is well resolved, it 
leads to the requirement that particle jitter be eliminated. 
To achieve this, the artificial viscosity must be calibrated 
to produce the correct results in test problems where the 
answer is known. Additionally, the widely discussed require¬ 
ment that artificial viscosity heating be small compared to 
the imposed cooling must also be satisfied. 

Unfortunately, satisfying all these requirements (3D 
simulations, no particle-jitter, minimal artificial viscosity 
heating) requires a very high particle number. Even for 
the parameters used in this paper, which were chosen to 
maximise computational efficiency, achieving an artificial 
viscosity heating less than 10% of the cooling rate when 
/3 = 20 with the standard fixed artificial viscosity method 
(ctspH = 1.0) requires of order 400 million particles in 3D. 

Despite these severe numerical limitations, the time 
scale arguments presented in this paper allow us to pre¬ 
dict that Pcrit can be at most a few tens. Given the strong 
scaling of /3 with radius in realistic protoplanetary discs, a 
more precise constraint on /3 C Ht is unlikely to be of physical 
importance for disc fragmentation as a mechanism for giant 
planet formation. 


7 MATERIALS & METHODS 

In the interests of reproducibility and transparency, all 
code and data used in performing this work have been 
made freely available online at https://bitbucket.org/ 
constantAmateur/discfragmentation 


8 ACKNOWLEDGEMENTS 

This work benefited greatly from discussions with Giuseppe 
Lodato, Deborah Sijacki, Richard Nelson, Farzana Meru, 
Richard Booth & Daniel Price on issues of both physics 
and numerics. We would like to thank Richard Booth & 
Sijme-Jan Paardekooper for providing valuable feedback on 
the manuscript. Sijme-Jan Paardekooper deserves special 
thanks for having provided advice and feedback throughout 
the course of this project. We thank the referee for comments 
which improved the paper. 

Matthew Young gratefully acknowledges the support of 
a Poynton Cambridge Australia Scholarship. The simula¬ 
tions in this work were performed using HPC resources al¬ 
located through DiRAC project DP022. This work has been 
supported by the DISCSIM project, grant agreement 341137 
funded by the European Research Council under ERC-2013- 
ADG. 


REFERENCES 

Binney J., Tremaine S., 2008, Galactic Dynamics: Second 
Edition. Princeton University Press 
Boss A. P., 2000, ApJ, 536, L101 
Clarke C. J., 2009, MNRAS, 396, 1066 
Clarke C. J., Harper-Clark E., Lodato G., 2007, MNRAS, 
381, 1543 

Cossins P., Lodato G., Clarke C., 2010, MNRAS, 401, 2587 
Cossins P., Lodato G., Clarke C. J., 2009, MNRAS, 393, 
1157 

Cullen L., Dehnen W., 2010, MNRAS, 408, 669 
Dehnen W., Aly H., 2012, MNRAS, 425, 1068 
Fabrycky D. C., Murray-Clay R. A., 2010, ApJ, 710, 1408 
Forgan D., Rice K., Cossins P., Lodato G., 2011, MNRAS, 
410, 994 

Gammie C. F., 2001, ApJ, 553, 174 
Hopkins P. F., Christiansen J. L., 2013, ApJ, 776, 48 
Kratter K. M., Murray-Clay R. A., 2011, ApJ, 740, 1 
Lattanzio J., Monaghan J., Pongracic H., Schwarz M., 
1986, SIAM Journal on Scientific and Statistical Com¬ 
puting, 7, 591 

Lodato G., 2008, New Astron.Rev., 52, 21 
Lodato G., Clarke C. J., 2011, MNRAS, 413, 2735 
Lombardi J. C., Sills A., Rasio F. A., Shapiro S. L., 1999, 
Journal of Computational Physics, 152, 687 
Marois C., Macintosh B., Barman T., Zuckerman B., Song 
I., Patience J., Lafreniere D., Doyon R., 2008, Science, 
322, 1348 

Masset F., 2000, A&AS, 141, 165 
Meru F., Bate M. R., 2011, MNRAS, 411, LI 
Meru F., Bate M. R., 2012, ArXiv e-prints 
Monaghan J. J., 1997, Journal of Computational Physics, 
136, 298 

Muller T. W. A., Kley W., Meru F., 2012, A&A, 541, A123 
Paardekooper S.-J., 2012, MNRAS, 421, 3286 
Paardekooper S.-J., Baruteau C., Meru F., 2011, MNRAS, 
416, L65 

Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347 
Rice W. K. M., Armitage P. J., Bate M. R., Bonncll I. A., 
2003, MNRAS, 339, 1025 

Rice W. K. M., Forgan D. H., Armitage P. J., 2012, MN¬ 
RAS, 420, 1640 

Rice W. K. M., Lodato G., Armitage P. J., 2005, MNRAS, 
364, L56 

Rice W. K. M., Paardekooper S.-J., Forgan D. H., Armitage 
P. J., 2014, MNRAS, 438, 1593 
Springel V., 2005, MNRAS, 364, 1105 
Thacker R. J., Tittley E. R., Pearce F. R., Couchman 
H. M. P., Thomas P. A., 2000, MNRAS, 319, 619 
Toomre A., 1964, ApJ, 139, 1217 


© 0000 RAS, MNRAS 000, 000-000 


